Comprehensive analysis of the coding and non-coding RNA transcriptome expression profiles of hippocampus tissue in tx-J animal model of Wilson's disease

Wilson's disease (WD) is an autosomal recessive disorder with a genetic basis. The predominant non-motor symptom of WD is cognitive dysfunction, although the specific genetic regulatory mechanism remains unclear. Tx-J mice, with an 82% sequence homology of the ATP7B gene to the human gene, are considered the most suitable model for WD. This study employs deep sequencing to investigate the differences in RNA transcript profiles, both coding and non-coding, as well as the functional characteristics of the regulatory network involved in WD cognitive impairment. The cognitive function of tx-J mice was evaluated using the Water Maze Test (WMT). Long non-coding RNA (lncRNA), circular RNA (circRNA), and messenger RNA (mRNA) profiles were analyzed in the hippocampal tissue of tx-J mice to identify differentially expressed RNAs (DE-RNAs). Subsequently, the DE-RNAs were used to construct protein–protein interaction (PPI) networks, as well as DE-circRNAs and lncRNAs-associated competing endogenous RNA (ceRNA) expression networks, and coding-noncoding co-expression (CNC) networks. To elucidate their biological functions and pathways, the PPI and ceRNA networks were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. A total of 361 differentially expressed mRNAs (DE-mRNAs), comprising 193 up-regulated and 168 down-regulated mRNAs, 2627 differentially expressed long non-coding RNAs (DE-lncRNAs), consisting of 1270 up-regulated and 1357 down-regulated lncRNAs, and 99 differentially expressed circular RNAs (DE-circRNAs), consisting of 68 up-regulated and 31 down-regulated circRNAs, were observed in the tx-J mice group when compared to the control mice group. Gene Ontology (GO) and pathway analyses revealed that DE-mRNAs were enriched in cellular processes, calcium signaling pathways, and mRNA surveillance pathways. In contrast, the DE-circRNAs-associated competing endogenous RNA (ceRNA) network was enriched for covalent chromatin modification, histone modification, and axon guidance, whereas the DE-lncRNAs-associated ceRNA network was enriched for dendritic spine, regulation of cell morphogenesis involved in differentiation, and mRNA surveillance pathway. The study presented the expression profiles of lncRNA, circRNA, and mRNA in the hippocampal tissue of tx-J mice. Furthermore, the study constructed PPI, ceRNA, and CNC expression networks. The findings are significant in comprehending the function of regulatory genes in WD associated with cognitive impairment. These results also offer valuable information for the diagnosis and treatment of WD.

Genome-wide identification of RNAs in mice hippocampus. We analyzed DE RNAs using the EdgeR method, following the criteria P < 0.05 and |log2FC|> 1.5. Venn diagram, volcano plot, and clustering map were used to show the DE RNAs, as shown in Fig. 3. Figure 3a-c indicate the Venn diagram, volcano plot, and clustering heatmap of DEGs.   Validation of the DECs using qRT-PCR. Since ninety-nine circRNAs were found to have significant changes with WD in RNA-seq, the reliability of the results is unknown. To confirm the reliability, six circRNAs with relatively high-level fold changes in RNA-seq were randomly selected for qRT-PCR validation, including www.nature.com/scientificreports/ three up-regulated circRNAs (mmu_circ_0001700, mmu_circ_0001662, and mmu_circ_00013859) and three down-regulated circRNAs (mmu_circ_0001859, mmu_circ_0000626, and mmu_circ_0000460). As presented in Fig. 4, the qRT-PCR data show that the circRNAs mmu_circ_0001700 and mmu_circ_0001662 were significantly up-regulated, while the circRNA mmu_circ_0001859, mmu_circ_0000626, and mmu_circ_0000460 were significantly down-regulated in the WD group. These results were consistent with RNA-seq, suggesting that RNA-seq data are available.    www.nature.com/scientificreports/ (GO: 0005515), and haptoglobin binding (GO: 0031720) (Fig. 5b). KEGG pathways enriched in mRNA surveillance pathway (mmu03015), Calcium signaling pathway (mmu04020), African trypanosomiasis (mmu05143), and Dopaminergic synapse (mmu04728) et al. (Fig. 5c). Reactome enrichment was shown in Fig. 5d, such as  www.nature.com/scientificreports/ p75NTR recruits signaling complexes (R-MMU-209543), NF-kB is activated and signals survival (R-MMU-209560), and NRIF signals cell death from the nucleus (R-MMU-205043). All enrichment results were displayed in Supplementary Table S5. Two KEGGs related to hippocampal neurophysiology and pathology, with DEGs indicated, are shown in Fig. 6, including Long Term Potential (mmu04720) (Fig. 6a), and Glutamatergic Synapse (mmu04724) (Fig. 6b).
Construction of ncRNA-miRNA-mRNA ceRNA network. Competitive endogenous RNA (ceRNA) regulatory networks play a crucial role in various human diseases. We constructed a circRNA-miRNA-mRNA and a lncRNA-miRNA-mRNA triple regulatory network based on the above differentially expressed ncRNAs. The functional parts and potential mechanisms of this network were assessed by functional enrichment analysis.

Discussion
The current work is the first transcriptome analysis of circRNAs, lncRNAs, and mRNAs in the hippocampus of the tx-J mouse model of WD. In these models, differentially expressed circRNAs, lncRNAs and mRNAs are associated with the mechanisms by which cognitive impairment occurs in WD. We focused on the participating components of non-coding RNAs in the ceRNA network and their effects on cellular functions. We analyze the characteristic involvement of non-coding RNAs in ceRNA networks regulating cognitive impairment in WD, which contributes to exploring biomarkers for the diagnosis or prognosis of cognitive impairment in WD at the molecular level. A growing body of evidence suggests that RNA transcription and metabolism changes are critical to developing complex symptoms and diseases. Non-coding RNAs represent the majority of the organism's transcriptome, which is abundant in the central nervous system and intensely involved in the transcriptional regulation of RNAs, thus being associated with the mechanisms of many neurological symptoms and diseases 17 . Several studies have shown that coding and non-coding RNAs regulate transcription through binding miRNAs and thus affect downstream target gene expression, which is known as a competitive endogenous RNA mechanism 18 . This pattern of gene regulation can be used to explain the complexity of species or phenotypic 19 .   20 , including motor, cognitive, and behavioral deficits, and that persistent or progressive cognitive dysfunction is present, even in patients treated with "copper detoxification" and patients with or without brain copper accumulation 21 . As mentioned previously, tx-J mice are currently an ideal model for studying the transcriptome sequencing of WD and related phenotypes. In this study, we confirmed the presence of cognitive impairment in tx-J mice through a water maze test, including decreased spatial learning and memory abilities. Meanwhile, regarding morphological examination of hippocampal structures, tx-J model mice showed destruction of neuronal cells in the dentate gyrus region of the hippocampus. They significantly decreased the number of new neuronal cells. It indicates that there is detectable cognitive impairment in tx-J mice compared with controls. On the one hand, this evidence of cognitive impairment can improve animal models' scientificity and strengthen sequencing results' reliability.
In this study, we performed the differential expression analysis of mRNA\circRNA\lncRNA, constructed a PPI expression network, and the ncRNA-related ceRNA network and ncRNA/mRNA co-expression network in the hippocampal tissue of WD model tx-J mice for the first time. We screened 361 DEGs. Many of these significantly www.nature.com/scientificreports/ differentially expressed target genes are closely associated with cognitive impairment pathology. We found that some mRNAs that miRNA targets in ncRNA-associated ceRNA regulatory networks were also differentially expressed mRNAs, such as Fosb, Shank3, and Cacna1i, etc. Fosb is an activity-dependent transcription factor that accumulates and is retained during chronic cellular activity due to its long half-life 22 . Fosb plays a critical role in regulating hippocampal memory; Studies have confirmed that the Fosb family regulates many addiction-related targets of crucial target genes through histone modifications 23 . It has also been demonstrated that over-expression of Fosb in the hippocampus may affect learning and memory 24 . A study examining Fosb regulation of gene expression in AD model mice with cognitive dysfunction suggests that Fosb may suppress the expression of c-Fos, an early gene critical for plasticity and cognition, by binding to promoters and triggering histone deacetylation and that the long half-life of Fosb makes it a possible cause of the persistence of cognitive deficits 25 . As an excitatory postsynaptic scaffolding protein, SHANK3 acts on various postsynaptic density proteins to regulate postsynaptic neurotransmitter receptors and signaling molecules 26 . A study based on the prediction of single-gene variant levels in Autism Spectrum Disorder (ASD) showed that SHANK3 mRNA levels affect synaptic transmission in the hippocampus of mice, leading to possible long-term depression and long-term impairment of learning and memory 27 . It was found that the voltage-gated calcium (Cav) channel gene CACNA1I (Cav3.3) is considered a risk factor for schizophrenia and that variants in this gene lead to disruption of neuronal excitability and brain network activity, affecting processes such as transmitter release, sensation, memory, and sleep 28 .
Non-coding RNA is a class of RNA molecules that do not encode proteins and account for approximately 98-99% of the total RNA in mammalian genomes 29 . There is increasing evidence that non-coding RNAs play important regulatory roles by actively interacting with other molecules; for example, a non-coding RNA interacts with one or more target molecules to regulate cells or pathways to affect normal physiological or pathological processes, including the regulation of neurological diseases 30 . CircRNAs are a class of covalently closed, single-stranded cyclic molecules without 5'-caps and 3'-polyadenylated tails 31 , which are stable, abundant, and conserved 32 and have unique potential for medical research. In this study, we identified the circRNAs profile in hippocampus tissue with WD, and 99 significant DECs were screened out. In particular, mmu_circ_0001859 and mmu_circ_0000242 were strongly related to Axon guidance, Wnt signaling pathway, Calcium signaling pathway, MAPK signaling pathway, Focal adhesion, Neurotrophin signaling pathway, TGF-beta signaling pathway, Cell cycle, ErbB signaling pathway, Notch signaling pathway, p53 signaling pathway, Fc gamma R-mediated phagocytosis, etc. www.nature.com/scientificreports/ Studies have confirmed the biological functions of circRNAs, such as participating in transcriptional regulation in the nucleus, acting as endogenous adsorbers of miRNAs, templates for protein or peptide translation, and regulators of gene expression 33 . CircRNA is exceptionally abundant in the mammalian brain, which is upregulated at the overall level during neuronal differentiation, highly enriched in neurosynapses, and plays a vital role in the development of neuropsychiatric disorders 34 . A study based on deep RNA sequence analysis systematically elucidated the cranial circRNA-associated ceRNA mechanism in an AD model mouse (APP/PS1 mice) and found that the circRNA-associated ceRNA network in this model mouse is mainly involved in dendritic development and memory (Sorbs2) and mouse neurodevelopment (ALS2), which provides new ideas for the clinical diagnosis and treatment of AD 35 . The first exploration of hippocampal circRNA expression profiles in aged mice with POCD suggests that mmu_circRNA_22058 and circRNA_44122\Egfr ceRNA network or circRNA_22673\Prkacb ceRNA network may have meaningful involvement in the course of POCD disease 36 . In this study, we identified 99 significantly differentially expressed circRNAs in hippocampal tissue of WD model tx-J mice compared to controls. We constructed a ceRNA network based on differentially expressed circRNAs and performed functional enrichment analysis. We performed GO and KEGG analyses on both up-and downregulated genes of circRNA-associated ceRNA networks respectively, and identified a series of enriched terms associated with neurological diseases and cognitive processes. In our results, GO of up-regulated genes was enriched in the cytoplasm (GO: 0005737), protein binding (GO: 0005515), nucleus (GO: 0005634), and metal ion binding (GO: 0046872). In contrast, KEGG was enriched in the Ras signaling pathway (mmu04014), PI3K-Akt signaling pathway (mmu04151), Calcium signaling pathway (mmu04020), etc. GO of down-regulated genes were enriched in cell junction (GO: 0030054), cell projection (GO: 0042995), synapse (GO: 0045202), and cytoplasmic vesicle (GO: 0031410), and KEGG was enriched in Amphetamine addiction (mmu05031), Adrenergic signaling in cardiomyocytes (mmu04261), Tight junction (mmu04530), Circadian entrainment (mmu04713). We identified several ceRNA networks that may contribute to the progression of cognitive impairment in WD. For example, in the ceRNA network with mmu-miR-6931-5p as the core and Fosb as the target gene, the upstream 44 circRNAs may bind mmu-miR-6931-5p as a sponge to regulate the expression of Fosb target genes.
Unlike circRNAs, LncRNAs are langer than 200 nucleotides in length that have 5' end caps and 3' end multimeric tails so that they are easily recognized by oligonucleotide-based RNA sequencing 37 . The expression of lncRNAs in the brain is tissue-specific 38 , and studies have confirmed the presence of regional and cell-specific expression patterns of lncRNAs in cognitively relevant memory brain regions such as the hippocampus, prefrontal cortex, and amygdala 39 , such as lncRNA Gm9968, which is significantly enriched in mouse hippocampus tissue,   40 . A study on the mechanism of LncRNA MEG3 action on cognitive function in AD model rats showed that upregulation of LncRNA MEG3 inhibited neuronal damage, reduced Aβ positive expression, and improved inflammatory indexes in AD rats, improving cognitive function as well 41 . In mouse models of vascular cognitive impairment, LncRNA TUG1 can bind and interact with BDNF proteins, and its overexpression can lead to cognitive dysfunction in mice after VCI 42  The study presented in this paper has several limitations. Firstly, the small sample size may limit the generalizability of the findings. Additionally, results obtained from animal models may not be entirely representative of the transcriptional profile of all phenotypes of Wilson's disease (WD), including female mice and human patients. Secondly, the assessment of spatial learning and memory in mice may not have been comprehensive enough. Furthermore, this study solely investigated the transcriptional profile of the hippocampus in the WD animal model, and other copper deposition tissues were not examined, which means that it was not possible to investigate gene expression in different tissues. Moreover, the study's results lack experimental validation, including www.nature.com/scientificreports/ genes with significant expression differences and predictive network models, which requires robust experiments to provide evidence for the results. Although our study sheds light on the mechanisms of cognitive impairment in WD based on transcriptomics, we acknowledge that the expression regulation mechanism of transcription molecules is complex and requires further research, such as the regulation of epigenetic modification.
In conclusion, investigating the pathological phenotypic mechanisms of diseases is a challenging and timeconsuming undertaking. However, the advancement of high-throughput sequencing technology has significantly facilitated the exploration of molecular mechanisms underlying diseases. Our examination of the protein-protein interaction (PPI) network and non-coding RNA competing endogenous RNA (ncRNA-ceRNA) network in the hippocampal tissue of tx-J mice, a model of Wilson's disease (WD), carries important implications for identifying potential biomarkers and target values essential for diagnosing, treating, and developing drugs for cognitive impairment in WD.

Methods
Animals. The study utilized ten male tx-J mice (C3HeB/FeJ-Atp7b tx-J/J) weighing 20 ± 2 g and ten male wild-type mice (C3HeB/FeJ) aged 8-10 weeks, acquired from the Jackson Laboratory through Beijing Vital River Laboratory Animal Technology, Ltd. All mice were individually housed in cages with controlled humidity (50-70%), room temperature (18-22 °C), and separate air supplies, with ad libitum access to food and water under an alternating 12 h light/dark cycle for eight weeks. For this study, there was an overlap of animal subgroups; three of the four tx-J mice underwent high-throughput sequencing after MWM testing, four of the remaining seven mice underwent hippocampal histopathology testing, and the remaining three underwent qRT-PCR experiments. The wild-type mice were also treated according to this protocol. Estrogen has been shown to influence neuroplasticity in several brain regions, regulate and mediate neuronal synapses and hippocampal formation, and affect learning and memory 43 . Thus, male mice were chosen as research subjects to minimize these effects.
The methods used in this study were performed in accordance with the relevant guidelines and regulations of the ARRIVE2.0 guidelines 44  Morris water maze (MWM) test. The Morris water maze experiment is routinely used to assess spatial learning and memory abilities in rodents, which consisted of a black circular pool (50 cm high and 90 cm in diameter) filled with 45 cm high water at 21-23 °C and mixed with a white non-toxic, odorless dye, in which a black circular table fixed 1 cm underwater in the center of the fourth quadrant. Four mice from each of the two groups were selected for the experiment. In the positioning navigation experiment, the two groups of mice were put into the water from the 45°angle position facing the pool wall in any of the four east, west, south, and north quadrants and tested four times a day. If the mice climbed up to the hidden platform within the 60 s and stayed on the platform for more than 5 s, they were considered to have found the platform. The time from water entry to platform finding was recorded as the escape latency. If the mice did not find the platform within the 60 s, they were guided to the platform and stayed there for 10 s, and the escape latency was recorded as 60 s. Compare the escape latency of the opposite platform quadrant, the average escape latency of three quadrants except for the platform quadrant, the swimming distance of the opposite platform quadrant, and the total swimming distance except for the platform quadrant on the fifth day between the two groups.
Hippocampus histopathology. Four mice from each group were selected for hippocampus histopathology. Mice were anesthetized by intraperitoneal injection of sodium pentobarbital (2 mL/kg; Shanghai Chemical Reagent Company) after fasting for 12 h and taken bilateral hippocampus tissues, fixed with 4% paraformaldehyde, dehydrated in ethanol and xylene, embedded in clear paraffin, sectioned, delayed, stained with hematoxylin-eosin, further dehydrated and sealed, and observed under a light microscope.
RNA extraction and high-throughput RNA sequencing. After isolating hippocampus tissue in the WD and control groups on ice boxes, we performed deep-sequencing of ribosome RNA samples from the hippocampus of three control and three tx-J mice.
Total RNA of Peripheral Blood Mononuclear Cells (PBMC) was prepared using MiRNeasy Mini Kit (Qiagen, Germany). Afterward, the RNA Clean XP Kit (Cat#A63987, Beckman Coulter, USA) and RNase-Free DNase Set (Cat#79254, Qiagen, Germany) were used to purify the total RNA. The purified RNA was quantified using two instruments for detecting RNA content and quality: NanoDrop 2000 (Thermo Fisher Scientific, USA) and an Agilent Bioanalyzer 2100 (Agilent Technologies, USA). According to the manufacturer's instructions, the libraries were prepared using the truseq® chain total RNA sample preparation kit (Illumina, USA). Qubit 2.0 fluorometer was used to quantify the library. Furthermore, these libraries were verified by Agilent 2100 Bioanalyzer. After confirming the size and molar concentration of the inserted fragment, the libraries were diluted to 10 pm with CBOT to generate clusters and sequenced on Illumina HiSeq 2500 (Illumina, USA). The library was constructed and sequenced in OG Biotech Inc (Aoji Biotech, Shanghai, China).
Identification and qualification of the expression level of RNAs. By using FastQC (v. 0.11.3, http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc) 45 , the RNA-Seq reads were quantified. Seqtk (https:// github. com/ lh3/ seqtk) removes the adapter sequence of Illumina TruSeq, mm10 ribosomal RNA reads, and lowquality reads. Then, using Hisat2 (version: 2.0.4), trimmed reads were mapped to the mouse reference genome (mm10) downloaded from Ensembl. The count of each gene in aligned reads was detected by StringTie (version 1.3.0). With Perl script, normalization of gene counts and fragments per kilobase of transcript per million mapped reads (FPKM) was performed. www.nature.com/scientificreports/ For circRNAs, All valid sequencing data were processed using the following steps: 1) The data were aligned to the mouse genome (http:// genome. ucsc. edu/) 46 by BWA-MEM software (version: 2.0.4); 2) The reads that were contiguously aligned to the genomes were discarded; 3) The unmapped reads were analyzed explicitly for backsplice junction sites using CIRI algorithms to identify the possible circRNAs 47 and counts were normalized by SRPBM (Spliced Reads Per Billion Mapping). Genes generating individual circRNAs were identified by matching the genomic locations of circRNAs to those detected by TopHat/Cufflinks using BED tools 48 . The conservation of circRNA between humans and mice was analyzed using the UCSC Liftover tool. It was utilized to map the 3´ and 5´ flank coordinates of each dysregulated circRNAs to the human genome coordinates. The criterion for homology was detecting the splice sites within ± 2 nucleotide distance around the putative human sites 49 .
For the identification of lncRNAs, single lane paired reads of 2 × 150 bp (PE150) in length were assembled using the transcript assembly software StringTie to remove known mRNAs and transcripts < 200 bp in length. The prediction software: Coding Potential Calculator (https:// github. com/ bioco der/ cpc) and Coding Non-Coding Index (https:// github. com /www-bioinfo-org/CNCI#install-CNCI) for lncRNA prediction of the remaining transcripts were then filtered from the lncRNA dataset. The expression levels of lncRNAs were measured by fragments per kilobase of exon model per million mapped reads (FPKM), and the expression abundance of known genes in different samples was calculated from FPKM values.
PPI network construction and analysis. The STRING database (https:// cn. string-db. org/) is an integrated database that affords information on interactions of protein-protein, not only including experimental interactions but also predicted interactions, which was used to predict the interaction of DE mRNAs with proteins. Cytoscape (v3.9.1) is used for constructing the PPI network. Using the MCAO plugin with the standard sets, the modules of the PPI network were identified with Cytoscape. GO, and KEGG pathway enrichment analysis was performed using the OECloud tools (https:// cloud. oebio tech. cn). Visualize DEGs and related paths through the Pathview website (https:// pathv iew. uncc. edu/). Construction and analysis of ceRNA. All sequences of mice-derived miRNAs were obtained from miR-Base version 21 52 . OG-Biotech's custom-built software, which is based on miRanda software and RNAhybrid, was adopted to predict miRNA targets that can interact with the DECs and DELs.
The number of predicted binding miRNAs was counted for all circRNAs. The mRNA target for predicted binding miRNAs was collected from miRdana 53 , targetScan 54 , and mirTarbase 55 database. The Venny was employed to obtain a strongly-supported mRNA target predicted by at least two databases. DECs and mRNAs that share the same miRNA binding site represented circRNA-miRNA-mRNA interactions.
For the cis-acting role of lncRNAs acting on neighboring target genes, coding genes 10 kb/100 kb upstream and downstream of lncRNA were searched for, and their functions analyzed. For the trans-acting role of lncRNA, identified based on expression level, the expressed correlation between lncRNAs and coding genes was calculated with the R function "cor. test". The ceRNA network was visually displayed by Cytoscape software.
The hub genes are calculated and displayed using the relevant algorithms in the Cytohubba plugin 56  Gene ontology (GO) enrichment and KEGG pathway analysis. For GO analysis, the biological process (BP), cellular component (CC), and molecular function (MF) analyses were performed. To further understand the function of the DE circ-/lnc-/mRNAs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched were analyzed using OECloud tools (https:// cloud. oebio tech. com) to reveal the high-level function and utility of the transcriptome of tx-J mice, the such as cells, organisms, and ecosystems. Significant enrichment terms for GO and pathways were picked out according to the threshold of P < 0.05.
Construction and analysis of circRNA-/lncRNA-/mRNA co-expression network. Based on the matrix data of DELs/DECs/DEGs, the Pearson correlation coefficients between DEL/mRNAs and DEC/mRNAs were calculated, respectively, and the threshold value of |r|> 0.7 and P < 0.05 to obtain an up-/down-regulated coexpression list between DEL/mRNAs and DEC/mRNAs. The co-expression network visualization was calculated